%% 
% 'a' needs to be Nx2 array, where first column is temperature, second is
% growthrate
load('gr_data_all_20230413_v2.mat');
close all
a = gr_data_all;
T = unique(a(:,1));

G = [];
Gerr = [];
for k=1:6
    g = reshape(a(:,k+1),[length(a)/length(T) length(T)]);
    gmean = nanmean(g,1)';
    gstd = nanstd(g,1)';

    G = [G gmean];
    Gerr = [Gerr gstd];

end


%% Plot data
T1plot = 16;
T2plot = 50;

clrs = distinguishable_colors(6);
clrs = cbrewer('qual', 'Paired',6);
strainnames = {'MG1655', 'BW25113', 'CS109', 'BL21', 'dnaK', 'uup'};

ms = 6;
figure;

toplot = [1:4];

for k=toplot
    g_plot = G(:,k);
    gerr_plot = Gerr(:,k);
    errorbar(T,g_plot,gerr_plot,'ko', 'markersize', ms, 'markerfacecolor', clrs(k,:));
    hold on;

    %% Perform analyses and print fitting results
    Tfit = T(4:11);
    gfit = g_plot(4:11);
    gfit_err = gerr_plot(4:11);
    flm_linear = fitlm(Tfit, gfit, 'weights', gfit_err);
    plot(Tfit, flm_linear.feval(Tfit), 'color', clrs(k,:),'linewidth', 1)

    Arrhenius_fit = @(b,x) b(1).*exp(-b(2)./(x+273.15));
    b0 = [1E12 8000];
    flm_Arrhenius = fitnlm(Tfit, gfit, Arrhenius_fit, b0, 'Weights', gfit_err);
    plot(Tfit, flm_Arrhenius.feval(Tfit), 'color', clrs(k,:), 'linewidth', 1);

    disp(strainnames{k})
    flm_linear
    flm_Arrhenius
    
end

set(gcf, "Position", [0 0 400 300]);
set(gca, 'FontSize', 20);
xlabel('Temperature (°C)');
ylabel('Growth rate (1/h)');
box off;
legend(strainnames{toplot})
xlim([T1plot T2plot])